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0\ ; Abstract 

We study the role of demographic fluctuations in typical endemics as exemplified by the stochastic 

' r-| ' 
Q ■ SIRS model. The birth-death master equation of the model is simulated using exact numerics and 

(p : 

analysed within the linear noise approximation. The endemic fixed point is unstable to internal 
demographic noise, and leads to sustained oscillations. This is ensured when the eigenvalues (A) of 
^ • the linearised drift matrix are complex, which in turn, is possible only if detailed balance is violated. 

In the oscillatory state, the phases decorrelate asymptotically, distinguishing such oscillations from 
— ^ , those produced by external periodic forcing. These so-called quasicycles are of sufficient strength 

Q ' to be detected reliably only when the ratio |/m(A)/i?e(A)| is of order unity. The coherence or 
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regularity of these oscillations show a maximum as a function of population size, an effect known 

►^ , variously as stochastic coherence or coherence resonance. We find that stochastic coherence can be 

00 ' 

\^ ' simply understood as resulting from a non-monotonic variation of |/?Ti(A)/i?e(A)| with population 

O ■ 

fT^ ■ size. Thus, within the linear noise approximation, stochastic coherence can be predicted from a 
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purely deterministic analysis. The non- normality of the linearised drift matrix, associated with the 
violation of detailed balance, leads to enhanced fiuctuations in the population amplitudes. 
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I. INTRODUCTION 

Two and a half centuries ago, D. Bernoulli [1] used a nonlinear ordinary differential 
equation to study the effect of cow-pox inoculation on the spread of smallpox. This was 
one of the earliest examples of the mathematical study of epidemics. This field of study 
continues to hold the interest of the scientific community especially in the light of recent 
outbreaks of viral pandemics like SARS and HlNl. Kermack and McKendrick in their 
seminal paper [2] put forward the classic Susceptible- Infected-Recovered (SIR) model of the 
spread of epidemics which, like most early epidemic models, assumes a homogeneously mixed 
population. More recent work focuses on the geotemporal spread of epidemics, especially 
on model networks [3, 4]. However, homogeneous mixing models still prove to be useful 
[5] and have been used to study various outbreaks of diverse size, fatality and chronology. 
Examples range from the study of the plague in the village of Eyam in 1665-66 [6] to the 
Bombay plague of 1905-06 [2] and the influenza epidemic in an English boarding school [7]. 

Mathematical models like the SIR model are usually analysed deterministically and are 
only exactly valid when the size of the population under consideration is exceedingly large. 
Fluctuations due to finite population sizes or due to external causes can give rise to phenom- 
ena which cannot be captured by deterministic mean-field models and necessitates the use of 
stochastic models. Bartlett [8, 9] was one of the first to realise that a stochastic description 
was necessary to explain the periodic recurrence of measles, a phenomenon which could not 
be explained by deterministic models [10, 11]. Bartlett formulated [8] a stochastic version 
of the SIR model to describe the periodic recurrence of measles. 

The mechanism for the generation of sustained oscillations in population dynamics has 
been analysed within the stochastic framework [12] which concentrates on external fiuctua- 
tions as the noise source. However, finite-sized populations give rise to fiuctuations whose 
relative amplitude is of the order of the inverse of the square root of the size of the popu- 
lation [13]. The role played by this internal noise, arising out of demographic stochasticity, 
in the generation of sustained oscillations has been studied in a prey-predator model using 
a master equation approach by McKane and Newman [14]. They have used the expansion 
method due to van Kampen [15] in their analysis, which provides a systematic way of de- 
riving the phenomonological equations due to Bartlett [8]. Alonso et al. [16] used similar 
techniques in an open model of infectious diseases within the homogeneous mixing assump- 



tion, while Rozhnova and Nunes [17] applied systematic expansion to a closed epidemic 
model on networks, using a pair approximation. The oscillations generated and sustained 
by internal noise are called endogenous resonant quasicycles and are qualitatively different 
from stochastic oscillations forced by external periodicities which are exogenous [18]. The 
quality or coherence of these oscillations are intuitively expected to vary monotonically with 
the size of the population or equivalently, the relative noise amplitude. However, it has 
been observed in various theoretical models including the Fitz Hugh-Nagumo [19] and gene 
circuit models [20] that the regularity or coherence of oscillations is small for low and high 
noise amplitudes and reaches a maximum for an intermediate value. This phenomenon is 
called stochastic coherence or coherence resonance and has also been observed in optical 
laser experiments [21]. 

In this work we analyse the generation of quasicycles due to internal noise, as well as 
the non-trivial variation of the quality of oscillation with respect to population size, in a 
closed epidemic model under the homogeneous mixing assumption. The closed system is 
relevant in many epidemiological situations, for instance in boarding houses [7], or island 
communities, where no inflows or outfluxes occur. Further, the conservation of populations, 
as implied by a closed system, allows one to deal with a lower- dimensional problem. We 
exploit this in a systematic manner and show how the master equation can be marginalised 
using the conservation constraint. The existence of an endemic fixed point allows a two- 
stage linearisation procedure to be carried out on the model. The linear noise approximation, 
followed by a further linearisation about the endemic fixed point, reduces the model to the 
standard multivariate Ornstein-Uhlenbeck (OU) form. Exploiting the linear and Gaussian 
character of the multivariate OU process then allows for stochastic behaviour to be predicted 
from the deterministic part of the dynamics, in a spirit similar to the Onsager regression 
method of equilibrium statistical mechanics. 

Below we review (section II) the deterministic analysis of the SIRS model, emphasising 
the behaviour around the stable endemic fixed point. Perturbations about this fixed point 
decay either monotonically or in a damped oscillatory fashion. Stochastic analysis (sections 
III and IV), however, shows that demographic noise destabilises this endemic fixed point, 
generating and sustaining oscillations. We show the existence of stochastic coherence (section 
V) by analytical means and then confirm it numerically. We show from purely deterministic 
analysis that there is stochastic coherence if the absolute value of the ratio of the imaginary 



and real parts of the eigenvalues of the linearised drift matrix shows a maximum when 
scanned against noise amplitude. The position of this maximum gives the population size 
corresponding to stochastic coherence for the relevant parameter values. We also show that 
(section VI) it is not possible to observe endogenous quasicycles unless detailed balance is 
violated. Finally we look at the non-normal aspect of the governing dynamics (section VII) 
and observe that the fluctuation amplitudes of the populations increase due to non-normality. 

II. SIRS LINEAR DETERMINISTIC ANALYSIS 

The classic SIR model for infectious diseases (S stands for susceptibles, I for infected and 
R for recovered) considers the population to be homogeneously mixed and constant in total 
number [2]. The SIRS model is a variant of the SIR model where the recovered section of the 
population lose their immunity after a delay and become susceptible. The nonlinear ODE 
system of the form n = f (n), where n = {S, /, R}, describing the SIRS model is constrained 
by the fixed population size Q and is hence a closed system. 

S = aR-(3SI 

i = ^SI-^I (1) 

R = ^I-aR 

The rate of infection is /3, the rate of recovery is 7 while a is the rate of loss of immunity. 
The fixed point (n = n*) is given by 
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(3 (3 \ a + 'J J /3\a + 7 
The steady state with zero infected is not of interest in the present study. The fixed point 
is endemic with non-zero infected in the steady state (/* > 0) when the condition f3Q > 7 
is satisfied. 

Since there is a constraint in the system, S + I + R = Q, the 3x3 system is effectively 
a 2 X 2 system with R = Q — S — I. 

S = a{Q-S-I)- PSI (3) 

/ = psi-^i 

The dynamics of small perturbations, 5n = {5S^ 5/}, about the fixed point are described 
by the linear ODE system 5ri = A ■ (5n. Here Aij = dfi/dnj\n=n* is the Jacobian matrix at 



the fixed point and is given by 
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Its eigenvalues are 
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the real parts of which are always negative for an endemic steady state since PQ > 7. Hence 
the endemic fixed point is always asymptotically stable. Perturbations about the fixed point 
decay monotonically if the eigenvalues are purely real and in an oscillatory fashion if they 
are complex. These correspond, respectively, to overdamped and underdamped decay. In 
Figure (1), we plot both time traces and phase portraits of S and I showing the underdamped 
and overdamped cases. Figure (2) is a state diagram of the model, showing the ratio | ^K/ 1 . 
The region of complex eigenvalues, corresponding to underdamped decay, is bounded by the 



contours labelled by 



I J-m(A) I 
I Re{X} I 
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III. SIRS LINEAR STOCHASTIC ANALYSIS 



Relative fluctuations about the deterministic expected values vary as the inverse of the 
square root of the number of interacting entities and thus become important when the enti- 
tites are few in number. Often one finds that this is indeed the case in biological systems [13]. 
Our present study concerns populations where fluctuations due to demographic stochasticity 
cannot be ignored and mean-field deterministic analysis fails to capture its non-trivial con- 
tributions. It then becomes necessary to employ stochastic methods to reliably understand 
the role of fluctuations. 

We begin by writing down the birth-death master equation (ME) of the SIRS model. Let 
the state of the system at any time t be given by the vector n = [rii, ^2, ^^3) where n^ is the 
number of individuals in each class (z = 1 for S, z = 2 for I and z = 3 for R). The general 
birth-death ME is then [22] 
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FIG. 1: (Color online) Underdamped and overdamped decay of perturbations. The top two plots 
show underdamped decay with parameter values /3 = 0.0021, a = 0.1,7 = 1.0 and population size 
il = 1000 (where the Jacobian has complex eigenvalues). The bottom two plots show overdamped 
decay with parameter values /3 = 0.0021, a = 5.0, 7 = 1.0 and population size fi = 1000 (where the 
Jacobian has real eigenvalues). The S vs I plot for underdamped decay shows a spiral while that 
for overdamped decay does not. The former is a stable spiral while the latter is a stable node. 
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+ Y, {^a (n - r„)P(n - r„, t) - t-(n)P(n, t)} (6) 



Here P(n, t) is the conditional probability for the system to be in the state n given some 
fixed initial state, t^ and t~ are the birth and death rate terms and r^ is the vector denoting 
the change in the number of entities in the a-th reaction. For the SIRS model we have 
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(3nin2; tj = 7723; tj = aris; t^ - ,3 - ^3 
(-1,+1,0); r2 = (0,-l,+l); r3 = (+1,0,-1) 

The variable N(t) = ni{t) + 7^2 (^) + n^lt) is constrained by the fixed population size Q. 
Incorporating this constraint within the ME allows us to work with a 2 x 2 system. The 
partial time derivative of P{N{t),t) vanishes if N{t) = Q for all t. We marginalise with 
respect to one of the variables (here we choose n^) taking the population size as parameter: 
P{ni,n2,t\Q) = J2n3^ni+n2+n3,nP{'n,t). Modifying the birth terms and the state change 
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FIG. 2: (Color online) Plot of the absolute value of the ratio of the imaginary and real parts of the 
eigenvalues A of the linearised Jacobian matrix A against the dimensionless parameters j3/"f and 
a/7 and population size = 1000. The outermost white contours, labeled by |Im(A)/Re(A)| = 
0, enclose the region where the eigenvalues are complex, which is a necessary condition for the 
existence of quasicycles. The imaginary parts of the eigenvalues are zero in the outermost two 
regions of the plot. The inner white contours, labeled "5Pii(a;)/9a; = 0" and "9^22 ('^)/9w = 0", 
enclose the regions for which the PSD shows a peak. This is a sufficient condition for the existence 
of quasicycles. The innermost black contour, denoting |Im(A)/Re(A)| = 1 and marked as such, 
encloses the region where the quasicycles are of sufficient strength to be reliably detected. This 
region is labeled "|Im(A)/Re(A)| > 1", which is a necessary and sufficient condition for the reliable 
detection of quasicycles. Each condition presented above is stricter than the previous, leading to 
a nesting of regions of parameter space as regards the existence and detection of quasicycles. 

vectors appropriately, i.e. replacing 723 by fi — ni — 7^2 and writing the r^ as 2 x 1 vectors, 
we get the marginalised ME. 
dP{ni,n2,t I f2) 



dt 



/3(ni + l)(n2 - l)P(ni + 1, n2 - 1, 1 1 f]) 

+ 7(n2 + l)P(ni,n2 + l,t|l^) 
+ a {Q — rii — n2 + l)P{ni — 1, n2, t \ Q) 

- {/3nin2 + 7^2 + a (fi - ni - 722)} -P(ni, ri2, t \ fi) (8) 

The transition probability for the infection step is non-linear and as such the ME is not 
solvable analytically. However, it is possible to simulate the ME using the Doob-Gillespie 
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FIG. 3: (Color online) Numerical simulation of susceptibles overlaid on deterministic underdamped 
decay. Parameters are (3 = 0.0021, a = 0.1,7 = 1.0 and population size il = 1000. There are noise- 
induced oscillations in the stochastic case which are not seen in the deterministic analysis. The 
time period of oscillations is approximately 20 in units of I/7 (simulations have been performed 
after non-dimensionalisation) . This corresponds well with the frequency seen (Figure 4) in the 
PSD analysis for the same set of parameter values. 

stochastic simulation algorithm (SSA) [23-25]. This generates an exact sampled trajectory 
of the jump stochastic process described by the ME. We non-dimensionalise time by work- 
ing in units of I/7. Figure (3) shows the numerical simulation of the susceptibles using 
the SSA, compared with a deterministic solution of the ODE system. The demographic 
fluctuations induce and sustain approximate cycles in the populations, a feature absent in 
the deterministic model. 

In the absence of exact solutions, we try to characterise these fluctuations within an 
approximation method due to van Kampen [15] which replaces the jump process with a 
stationary multivariate Ornstein-Uhlenbeck process. The Gaussian nature of the process 
can then be utilised to obtain analytical solutions for the fluctuation properties, while the 
linear nature of the process can be utilised to make connections between the deterministic 
and fluctuating dynamics. 

We expand the variables in the population size Vl (the large parameter of the approxima- 
tion method) so that the size of the jumps decreases as the population is increased, 

n = l]n + f]^/^x (9) 

where n is the mean value of n and x denotes the fluctuations around the mean. Assuming 



the fluctuations obey a diffusion process about tlie mean yields a Fokker-Planck equation 
(FPE) for the fluctuations, 



dtPi^,t) = -a,, [Ai(x)P(x,t)] + -d,dj [5,,(x)P(x,t)] 



(10) 



where repeated indices indicate summation, dt = d/dt and di = d/dxi. This is the linear 
noise approximation. The elements of the drift vector A(x) and the diffusion matrix B(x) 
are given, following the prescription in Gardiner [22], as 

2 

^t (11) 

5.,(x) = ^rXt+(x) 



a=l 



{i,j = 1,2 being the component indices). Linearising a second time about the endemic fixed 
point we get the FPE of a stationary multivariate Ornstein-Uhlenbeck process 



9tP(x,t) = -^ 



«j 



Ajd, {x,P(x,t)} - -Bi,did,P{^,t) 



(12) 



where Aij and Bij are the elements of the linearised drift and diffusion matrices. For the 
SIRS model, their values are (from Equation (11) after putting x = x*) 
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We note that this linearised drift matrix A is identical to the linearised Jacobian matrix 
(Equation 4) obtained from the deterministic analysis and hence the two matrices share the 
same spectrum. This allows us to predict, under the two-stage linearisation procedure, the 
existence of non-trivial stochastic phenomena like noise-induced quasicycles and stochastic 
coherence purely from a deterministic analysis of the spectral structure of the linearised 
Jacobian. We shall discuss this important point in greater detail in sections IV and V. This 
also allows us to use the terms "linearised drift matrix" and "linearised Jacobian matrix" 
interchangeably. 



The multivariate Ornstein-Uhlenbeck process has exact solutions for both stationary and 
transition probability densities. Both are multivariate Gaussians, fixed by the equal time 
covariance matrix Sj^ = {{xiXj)) and the matrix of time correlations Cij{T) = {{xi{t)xj{t + 
r))), where the double angular brackets denote the cumulant [15]. S can be obtained by 
solving the steady state Einstein relation [15, 22]. 



AS + SA' + B = 



(15) 



This has the form of a matrix Lyapunov equation, and can be solved using a method first 
proposed by Barnett and Storey [26] in the context of linear control systems. We note that 
Equation (15) can be written as the sum of a matrix and its transpose S + S"'" = where S 
is the anti-symmetric matrix AS + |B. We can solve for S in terms of A and B using the 
relation 



AS + SA^ = - (BA^ - AB) 



(16) 



which is obtained by eliminating S from the Einstein relation and using the definition of 
S. Since S is anti-symmetric, it is specified by a single parameter when it is of size 2x2. 
This parameter can be obtained directly from Equation (16), since both A and B are two- 
dimensional matrices and are known. For higher dimensions, matrix decompositions are 
convenient when solving for S. 

For the SIRS model (using Equations 13 and 14), we have 

-1 



-fifSn - 7)(a^ + 2^7 + 27^ + af3n) \ 
2/3(« + 7)(« + /3fi) J 
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Knowing S, A and B we can now write down the covariance matrix 

1. 



which for the SIRS model is 
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a^+7^+a(/3Q+7) 
a{a+/3n) 

-1 



B 



a(a+/3n)^ +7(0+7) (/3n-7) 
/3(a+7)2(a+/30) 



Having obtained the matrix S, the matrix of time correlations follows as 



C(r) = ((x(i)x(i + r))) 



^tA-< 



(17) 



(18) 



(19) 



(20) 



The stochastic SIRS model, in the linear noise approximation, is completely specified by S 
and C(r). In the next section we use quantities derived from these to examine the model 
for signatures of oscillatory behaviour. 
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IV. NOISE-INDUCED OSCILLATIONS: ENDOGENOUS QUASICYCLES 

The trace of the variation of the populations with time shown in Figure (3) is strongly 
suggestive of sustained oscillations. This can be verified quantitatively by measuring the 
power spectral density (PSD) of the population time series. A peak in the PSD indicates 
the presence of oscillations. The PSD matrix, in terms of the linearised drift and diffusion 
matrices for a multivariate Ornstein-Uhlenbeck process is 

P{uj) = {-iujI + A)-'^B{iujI + A^)-^ (21) 

where I is the identity matrix. The diagonal elements of this matrix give an estimate of the 
periodicity in the relevant variables (here S and I). The Pa for the SIRS PSD are 



T^ + uP 



p-(^) = 2^ ( .,4;:,2,. ) (22) 



where d = ^^J^fci, Pi = (a + 7)^ r^ = aM«' + 7^ + /3fi(/3^ - 7) + «(/3fi + 7)} /(« + 7)', 
q = a {a\a + 27) - 2-f\l3n - 7) + a(/3fi - 27)^} /(« + 7)2 and r = a\f3n - 7)^. 

In Figure (4) we plot the PSD for both S and I, comparing numerical simulation with 
Equation (22). A peak is clearly visible for parameters corresponding to underdamped 
dynamics. The peak disappears for overdamped dynamics as shown in the inset. The peak 
frequency (around Up = 0.3) corresponds to the period (T = 20) of the numerical time-trace 
(Figure 3). The excellent agreement between numerics and analytics provides a post-facto 
justification of the linear noise approximation for this problem. 

The PSD has peaks at real frequencies if and only if the extremum condition dPii{uj)/duj = 
has real roots. The regions of parameter space for which this occurs are bounded by 
contours labelled ''dPn{u)/duj = 0" and ''dP22{(^)/du = 0" in Figure (2). These are 
sufficient conditions for the existence of quasicycles. This approach has been used previously 
in the literature to detect quasicycles [14, 16, 17]. 

While Fourier analysis of a signal is a natural tool for studying oscillatory behaviour, 
a corresponding time-domain analysis must yield equivalent results. The time-correlation 
function forms the basis of a time-domain analysis, which for the multivariate Ornstein- 
Uhlenbeck process is given by Equation (20). The temporal variation of the time correlation 
is fixed entirely by the drift A which is the deterministic part of the dynamics, while its 
scale is set by S which involves the stochastic part of the dynamics through B. Defining a 
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FIG. 4: (Color online) Normalised power spectral density for S and I. There is excellent agreement 
between the analytically calculated (dashed blue for S, solid red for I) and numerically computed 
(blue dots for S, red dots for I) PSDs, thus justifying the linear noise approximation. The main 
graph is for parameter values /? = 0.0021, a = 0.1,7 = 1.0 at population size Q = 1000 which 
falls within the underdamped zone. The PSD peaks around frequency oj = 0.3, which corresponds 
approximately to the time period of the numerical signal as well as that of the ACF for the same 
set of parameter values (see Figures 3 and 5). The inset shows the PSD for the same size of the 
population at parameter values /? = 0.0012, a = 2.5,7 = 1-0 which falls within the overdamped 
zone and does not show any peak. 

normalised time correlation c(r) = C(r)S~^, we find that c(r) = e^^. This is of the form 
c{t) ~ exp[/2e(A)r] sin[/??i(A)r)]. This observation motivates the use of the ratio | ^\J \ 
to reliably detect quasicycles within the linear noise approximation, where A = eig(A). If 
the decay time scale, fixed by Re(A), is too short compared to the oscillatory time scale 
fixed by Im(A), the decay will dominate and oscillatory effects will not be discernible. This 
will be so even when the extremum condition has real roots. We thus propose a condition 
for clearly discernible quasicycles, namely | ^\J \ > 1. In Figure (2) we plot the contour 
I ^\J I = 1. The region I-^^t^I | > 1 is bounded on the right by this contour. As this is 
more stringent than the extremum condition dPii{uj)/duj = 0, it is entirely contained by the 
regions where the PSD has a peak. In Figure (5) we emphasise this point by comparing the 
ACF when the PSD has peaks at finite frequencies. When | ™;,, | is small the oscillations 



Re{X) 



I Im(X) I 



are barely discernible as seen from the rapid decay of the ACF. For I^^t^I of order unity 
clear signatures of oscillation are visible in the ACF. 
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FIG. 5: (Color online) Normalised autocorrelation of susceptibles and infected for parameter values 



P = 0.0021, a = 0.1, 7 = 1.0 and population size Q = 1000 which falls within the 



Re{X) I 



> 1 zone. 



The thick black line is the x-axis. There is clear oscillatory decay with a period of approximately 
20 in units of I/7 which agrees well with Figures (3) and (4) for the same set of parameter values. 
The inset plot shows the ACF for parameter values /3 = 0.009, a = 1.2, 7 = 1.0 and population size 
0, = 1000 which falls within the zone bounded by the contours dPii{u})/duj = 0, dP22{^) / duo = 



and 



I Im{X) I 



-p^T^l — 0, i.e the region where susceptibles should cycle according to the PSD analysis. The 
thick black line is once again the x-axis. There is a single zero-crossing of the ACF for S, which 
indicates non-oscillatory decay [18]. This shows the unreliability of PSD analysis in detecting 
quasicycles. These analytical plots have also been compared with numerical data (not shown here) 
with good agreement. 



V. QUALITY OF NOISE-INDUCED OSCILLATIONS: STOCHASTIC COHER- 
ENCE 



Noise-induced oscillations, unlike genuine oscillations, are not phase coherent and as such 
are called quasicycles. The coherence or regularity of quasicycles can be quantified by several 
measures. Here we use the quality factor, which measures the sharpness of the peak of the 
PSD, and the coefficient of variation which measures the regularity of the zero crossing of 
the signals. 

The quality factor, Q, is a dimensionless parameter that characterizes an oscillator's 
bandwidth relative to its peak frequency. 



Q = ujp/Au 



(23) 
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where Up is the peak frequency and Auj is the bandwidth. A high Q corresponds to oscilla- 
tions of greater regularity. We calculate the Q for the diagonal entries of the PSD matrix. 
Let ki be half the maximal power ki = ^Pu^ul) for each i = 1,2. We calculate the bandwidth 
or the full-width at half-maximum (FWHM) using the ki and Equation (22) to get 



(Aw), = yji2d/ki -q)- 2^r - 2dTi/ki (24) 

Using Equation (22), the peak frequencies are given by the positive square roots of the 
positive roots of the two quadratic equations z"^ + 2TiZ + (Fjg — r) = (for z = 1, 2) where 
z = iS^ and Fj, q and r are as defined in the previous section. The peak frequency and 
the FWHM together give the Q. Figure 6(a) shows a scan of the quality factor against 
population size and against the inverse of population size (inset). As one would expect, 
Q is low for high noise amplitudes and starts increasing as the noise decreases, keeping in 
mind that the relative noise amplitude varies as the inverse of the square root of the size 
of the population. However, the graph then has a maximum and then decreases for high 
amplitudes of noise. This is stochastic coherence. 

The coefficient of variation, Cy, is the variance over the mean of the times T between 
succesive zeros of a temporal signal. A sharp peak in the histogram of the intervals between 
zero crossings, then, indicates a strongly coherent signal. Cy is a dimensionless measure of 
this, 

C. . ^^^ (25) 

A low Cy indicates a high degree of coherence in the signal. Similar measures are used in 
the literature (for example [19] and [20]). Figure 6(b) shows the Cy for the mean crossing 
time of the numerical signal of the susceptibles scanned against population size and (inset) 
against its inverse . The plot shows a minimum which indicates stochastic coherence and 
hence numerically supports the analytical result given by the Q. 

Although this non-intuitive variation of the coherence of the quasicyles with the size of 
the population has a stochastic origin, it is controlled purely by the deterministic part of the 
dynamics. The analysis using the Q and the Cy require a knowledge of the diffusion matrix 
B. However, after the two-step linearisation procedure, the entire non-trivial dependance 
on the population size is contained only in the spectrum of the linearised drift matrix, while 
the diffusion matrix scales linearly with fi, as given by Equations (13) and (14). Thus, 
any non-monotonicities in the fluctuations arise purely from the deterministic part of the 
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FIG. 6: (Color online) (a) Quality factor and (b) Coefficient of Variation for the susceptibles 
(plotted for 200 runs) against population size (fi) and (inset) against inverse of population size. The 
solid line in (b) is a fifth-order polynomial fit. Parameter values are /3 = 0.0021, a = 0.1, 7 = 1.0 for 
(a) and j3 = 0.002, a = 0.1, 7 = 1.0 for (b). In the Q plot, the peaks for S and / are respectively at 
il = 1040 and il = 1000, while the Cv peak for 5, computed from the minimum of the fifth-order 
polynomial fit, is at = 985. The Cy plot for the infected (not given here) shows a peak (in the 
sixth-order fit used there) at $7 = 1135. 

dynamics, while the noise merely excites these modes. For a system which can be reduced to 
a standard multivariate Ornstein-Uhlenbeck process, the linearised drift matrix is identical 



to the linearised Jacobian matrix. This motivates the use of the ratio 



I Im{X) I 
I Re(X) I 



in determining 



the size of the population at which stochastic coherence is observed. This allows us to study 
stochastic coherence from the deterministic part of the dynamics. 
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FIG. 7: (Color online) Absolute value of the ratio of the imaginary and real parts of the eigenvalues 
of the linearised Jacobian matrix plotted against population size (Q) and (inset) against inverse of 
population size. Parameter values are /3 = 0.0021, a = 0.1,7 = 1.0. There is a peak at = 1000 
corresponding to the stochastic coherence point. 

We observe that in Figure (7) the ratio | ^^J \ , when scanned against the size of the 
population (Q), shows a peak which occurs at ^2 = 1000 for the parameters (3 = 0.0021, 
a = 0.1,7 = 1.0. We see that this value matches well with the peaks in Figures 6(a) and 
6(b), within numerical errors. We have calculated the peak value of the ratio in terms of 
the model parameters. If fip is the population size at stochastic coherence, then 

a + 27 



i Ln 



/3 



(26) 



Since the ratio is always positive, there is stochastic coherence for all values of parameters 
for which quasicycles exist. 

VI. DETAILED BALANCE VIOLATION NECESSARY FOR QUASICYCLES 



Quasicyles and stochastic coherence are not possible unless detailed balance is violated 
in the master equation. Typically, variables characterising biological systems are even under 
time-reversal, x{—t) = x{t), and this is true for the variables of the SIRS model. If the 
system is in state ng at time to = (say) and is in state n at some later time t, then 
the joint probability of the forward transition (no(to) — ^ ^{t)) is -P(n, t;no,0) while that 
of the reverse (n(to) — > iio(t)) is P(no,t;n, 0) provided all time-reversal parities are even. 
Microscopic reversibility implies that at equilibrium the steady state forward and reverse 
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joint probabilities must be equal. This is the condition for detailed balance and for a Markov 
process can be written as 

P(n, t|no, 0)P,(no) = P(no, t|n, 0)P,(n) (27) 

where the subscript 's' denotes steady state. Expressing this condition macroscopically in 
terms of the correlation function, expanding in Taylor series, and keeping the first order 
terms one obtains [27] the Onsager relations 

AS = SA"^ (28) 

which is the macroscopic condition for equilibrium. This condition requires the drift matrix 
to be related by a similarity transformation to a symmetric matrix [22] and hence restricts its 
spectrum to the real axis. Since it is not possible to have quasicycles without having complex 
eigenvalues, the violation of detailed balance is a necessary condition for the existence of 
noise-induced oscillations. 

Recalling that S = AS + |B and using the symmetry properties of S and B we can 
write down the following expression for S [28] 

S = i(AS - SA^) (29) 

which then is a measure of the deviation from detailed balance. The SIRS S matrix (Equa- 
tion 17) can never be zero for any choice of parameters under the endemic condition PQ > 7. 
Thus the SIRS model always violates detailed balance and therefore allows for quasicycles 
for any choice of parameters. 

VII. NON-NORMALITY INCREASES VARIANCE 

We have already noted that the violation of detailed balance is necessary for quasicycles. 
Here we further note that detailed balance violation has another consequence, that of en- 
hancement of fluctuation amplitudes. With detailed balance the drift matrix A is similar 
to a symmetric matrix, and is therefore normal (AA = A-'" A). In the absence of detailed 
balance, the drift matrix is no longer symmetric, and in this case is also non-normal. 

As has been noted by loannou [29], the variance of a non-normal system driven by diag- 
onal white noise is larger than its normal counterpart. Consider two stationary multivariate 
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FIG. 8: (Color online) Enhanced variance of time-traces of susceptibles and infected for parameter 
values P = 0.0021, a = 0.1, 7 = 1.0 with population size il = 1000. The solid black lines correspond 
to n lb \m, where n is the mean. 

Ornstein-Uhlenbeck processes with drift and diffusion matrices (Ai,B) and (A2, B) where 
Ai is non- normal but shares the same eigenvalues as the normal A2. Then, Schur decompo- 
sitions of the two matrices gives Ai = U(D + T)U^ and A2 = UDU^ where U is unitary, D 
is diagonal matrix of eigenvalues and T is strictly upper triangular. Restricting the forcing 
to be diagonally correlated white noise (B = I), loannou shows that Tr(I]i) > Tr(S2), where 
Si and S2 are the respective covariance matrices and Tr(. . .) denotes the trace of a matrix. 
For a general B which is not necessarily diagonal. Si = Af ^ (S — |B) and S2 = — |A^^B. 
We have calculated the ratio of the traces of Si and S2. 



Tr(Si) _ ^ Tr (A^^S + 2kUTU^B) 
IV(S2) Tr(S2) 



(30) 



where A is the determinant of Ai. This expression is valid only when the spectrum of Ai 
is purely real. For the SIRS model, this ratio is greater than unity. 

Individual time-traces also show an increase in variance. Figure (8) shows time-traces of 
S and I where the fluctuations are seen to be higher than the expected standard deviation 
values (n ± y/^, where n is the mean) marked by the black lines. 
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VIII. CONCLUSION 

In this paper, we have analysed a closed endemic model for sustained, though asymp- 
totically incoherent, oscillations in the population classes. These oscillations are generated 
through fluctuations brought about by internal demographic stochasticity which destabilise 
the endemic fixed point. The closed nature of the problem allows one to deal with a sim- 
plified lower-dimensional problem, an aspect we have exploited systematically by showing 
how the master equation can be marginalised using the constraint. This model also lends 
itself to a two-stage linearisation procedure, at the end of which it is reduced to a multi- 
variate Ornstein-Uhlenbeck form. This results in the identification of the linearised drift 
matrix with the deterministic Jacobian matrix linearised about the endemic fixed point and 
permits the analysis of stochastic behaviour from the deterministic behaviour. 

Noise-induced oscillations or quasicycles are possible only if the eigenvalues of the lin- 
earised Jacobian matrix are complex. These oscillations are distinct from those produced 
by external periodic agencies because their phases decorrelate asymptotically. Quasicycles 
can be reliably detected only if the oscillation time period is at least of the same order as 
the decorrelation time scale, as otherwise the decay dominates over the oscillation. Strong 
quasicycles are seen when the imaginary parts of the eigenvalues are larger than the real 
parts. 

Stochastic coherence, or the non-trivial maximisation of regularity of the oscillations at 
intermediate relative noise amplitudes (or equivalently at intermediate population sizes) , is a 
striking aspect of the SIRS quasicycles. We have seen this both analytically from the relative 
strength of the peak of the power spectral density and numerically by directly computing 
the signal-to-noise ratio of the time-traces of each population class. This analysis requires a 
knowledge of the intrinsic noise in the system, namely the diffusion matrix B. However, we 
find that, for systems which can be reduced to a standard multivariate Ornstein-Uhlenbeck 
form by the two-stage linearisation procedure mentioned earlier, it is possible to predict 
stochastic coherence purely from the deterministic analysis. Any non-trivial dependance 
on population size is contained only in the eigenvalues of the linearised drift matrix or 
equivalently the linearised Jacobian matrix, while the diffusion matrix scales linearly with 
the population size. Thus, any non-monotonicities in fluctuations arise entirely from the 
deterministic part of the dynamics, i.e. the spectrum of the drift matrix. The noise merely 
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excites these modes. This motivates the maximisation of the ratio I J!^/J I in the investigation 
of the population size value at which stochastic coherence is observed. Numerical results 
support this procedure. Therefore, we conclude that it is possible to make predictions 
about non-trivial behaviour of such systems in the stochastic regime by simply analysing 
the linearised deterministic dynamics. 

The violation of detailed balance is a necessary condition for the existence of quasicycles. 
Analysis of the drift, diffusion and S matrices indicates that the population system described 
by the SIRS model is always out of equilibrium and allows for quasicycles about the endemic 
fixed point for any choice of model parameters. Violation of detailed balance due to the 
non-normal nature of the system dynamics is manifest in the enhancement of fluctuation 
amplitudes of the populations. We have given an expression for the ratio of the trace of the 
non-normal covariance matrix over its normal counterpart, restricted to parameter values 
where the Jacobian spectrum is purely real. Numerics indicate that this ratio is greater than 
unity for the SIRS model. 

The analysis of this paper shows that the phenomenon of noise-induced oscillations 
and stochastic coherence can generically be expected in non-equilibrium birth-death jump 
Markov processes which can be reduced to the standard multivariate Ornstein-Uhlenbeck 
form by a successive application of two linearisation procedures: the linear noise approxi- 
mation followed by a linearisation about the fixed point of the system. This may therefore 
explain the appearance of asymptotically incoherent oscilations in other systems described 
by such equations, as for instance, in the repressilator [30]. 
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